In [3]:
import datetime
station_list=[41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
#station_list=[42027]
station_list_cs=[43150, 42867, 43014, 42339, 40990, 40948]
date_min=datetime.datetime(2011,5,1,0,0,0)
date_max=datetime.datetime(2011,10,1,0,0,0)
date_title = '%s_to_%s' % (date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'))
#Indexes - pressures = 1, ptemps = 2, vptemps = 3, rh = 4, uwind = 5, vwind = 6, calculated geopotential heights = 7 . . . may change
# Single level indexes - dates_for_plotting_single = 0, pwater = 1, lclpress = 2, lfcpress = 3, lnbpress = 4, cape = 5,cin = 6
In [17]:
########################################
# Read station radiosonde files in and get info on sounding times, dates, and heights reached
#
# Created by: Peter Willetts
# Created on: 14/5/2014
#
########################################
#
# http://www1.ncdc.noaa.gov/pub/data/igra/readme.txt
###################################################
import glob
import re
import numpy as np
# Create list of monthly radiosonde files to be used as input
# Derived parameters .dat files
#rad_flist = glob.glob ('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/derived_parameters/*.dat')
rad_flist = glob.glob ('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/*.dat')
match_header = re.compile(r'#*20110[56789]')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
# Read station name file
station_metadata=[]
f = open(station_list,'r')
for line in f:
line = line.strip()
line=re.sub(r'([A-Z])\s([A-Z])', r'\1_\2',line)
line=re.sub(r'([A-Z])\s\s([A-Z])', r'\1_\2',line)
station_metadata.append(line.split())
f.close()
# Read station files searching for start of each indpendent sounding header
all_stations_soundings=[]
for i in rad_flist:
sounding=[]
station_soundings=[]
station_no=i[-9:-4]
match_filename = re.compile('%s' % station_no)
for sm in station_metadata:
if match_filename.search(str(sm)) is not None:
station=sm[2]
latitude=sm[3]
longitude=sm[4]
#GUAN_code=sm[6]
#print station
#print latitude
#print i
f = open(i,'r')
for line in f:
line = line.strip()
columns = line.split()
#print line
# Extract independent header variable
if match_header.search(line) is not None:
#print line
date=columns[0][1:14]
time=columns[0][16:20]
time_hour=columns[0][14:16]
#print time_hour
#print columns[0][1:17]
no_levels=columns[1]
sounding=(date, time, time_hour, no_levels)
station_soundings.append(sounding)
all_stations_soundings.append((station, latitude, longitude, station_soundings))
np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived', np.array(all_stations_soundings, dtype=object))
#np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_derived', np.array(all_stations_soundings, dtype=object))
In [794]:
%matplotlib inline
import matplotlib
import matplotlib.mlab as ml
import scipy.interpolate
# Need to grid irregular radiosonde data before plotting
def grid_data(q, param):
xi=matplotlib.dates.drange(min(q[0]), max(q[0]), datetime.timedelta(hours=24))
yi=np.linspace(min(q[1]), max(q[1]), 25)
zi = ml.griddata(matplotlib.dates.date2num(q[0]),q[1],param,xi, yi, interp='nn')
#zi = scipy.interpolate.griddata((matplotlib.dates.date2num(q[0]),q[1]), param, (xi[None,:],yi[:,None]), method='linear')
return xi,yi,zi
def grid_data_cs(pressure, distance, param):
xi=np.linspace(0, max(np.array), 25)
yi=np.linspace(min(q[1]), max(q[1]), 25)
#yi=[1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20,10]
zi = ml.griddata(distance, pressure,param,xi, yi, interp='nn')
#zi = scipy.interpolate.griddata((matplotlib.dates.date2num(q[0]),q[1]), param, (xi[None,:],yi[:,None]), method='linear')
return xi,yi,zi
station_list_search='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
# Read station name files searching for start of each indpendent sounding header
station_metadata=[]
f = open(station_list_search,'r')
for line in f:
line = line.strip()
line=re.sub(r'([A-Z])\s([A-Z])', r'\1_\2',line)
line=re.sub(r'([A-Z])\s\s([A-Z])', r'\1_\2',line)
station_metadata.append(line.split())
f.close()
def station_name_search(stat):
for line in station_metadata:
if "%s" % stat in line: st = line[2].lower().title().replace('_',' ')
return st
def station_info_search(stat):
for line in station_metadata:
if "%s" % stat in line:
st = line[2].lower().title().replace('_',' ')
lo = float(line[3])
la = float(line[4])
return st,la,lo
def plot_rad(xi,yi,zi, min_contour, max_contour):
clevs = np.linspace(min_contour, max_contour,256)
ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
plt.figure(figsize=(12,8))
cmap=plt.cm.jet
cont = plt.contourf(matplotlib.dates.num2date(xi),yi/100, zi, clevs, cmap=cmap, extend='both')
cbar = plt.colorbar(cont, orientation='horizontal', pad=0.05, extend='both', format = '$%d$')
#cbar.set_label('$W m^{-2}$')
cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
cbar.set_ticklabels(['${%d}$' % i for i in ticks])
plt.gca().invert_yaxis()
plt.ylabel('Pressure (hPa)')
return cont,cbar
def plot_lcl_etc(xl,yl, label, colour, linestyle):
plt.plot(xl,yl, color=colour, label=label, linestyle=linestyle)
def grid_data_cs(pressure, distance, param):
xi=np.linspace(0, max(distance), 25)
yi=np.linspace(min(pressure), max(pressure), 25)
zi = ml.griddata(distance,pressure,param,xi, yi, interp='nn')
#print z
#zi = scipy.interpolate.griddata((matplotlib.dates.date2num(q[0]),q[1]), param, (xi[None,:],yi[:,None]), method='linear')
return xi,yi,zi
def plot_rad_cs(xi,yi,zi, min_contour, max_contour):
clevs = np.linspace(min_contour, max_contour,256)
ticks = (np.arange(min_contour, max_contour+tick_interval,tick_interval))
plt.figure(figsize=(14,8))
cmap=plt.cm.jet
cont = plt.contourf(xi,yi/100, zi, clevs, cmap=cmap, extend='both')
cbar = plt.colorbar(cont, orientation='vertical', pad=0.05, extend='both', format = '$%d$')
cbar.set_label('$W m^{-2}$')
cbar.set_ticks(np.arange(min_contour, max_contour+tick_interval,tick_interval))
cbar.set_ticklabels(['${%d}$' % i for i in ticks])
plt.gca().invert_yaxis()
plt.ylabel('Pressure (hPa)')
plt.xlabel('km from first station')
return cont,cbar
def calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon):
fslat_rad = radians(first_station_lat)
fslon_rad = radians(first_station_lon)
lat_rad = radians(station_lat)
lon_rad = radians(station_lon)
#Haversine Formula
a = sin((lat_rad-fslat_rad)/2)**2 + cos(lat_rad) * cos(fslat_rad) * sin((lon_rad-fslon_rad)/2)**2
c = 2 * atan2(sqrt(a), sqrt(1-a))
d = 6371 * c
return d
In [39]:
#################
# Plot frequency of radiosonde soundings by station vs date/time
##########
import os
import cPickle as pickle
import datetime
import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
from collections import defaultdict
import numpy as np
import matplotlib.cm as cm
# Load file from radiosonde_read_separate.py
import re
from matplotlib import rc
from matplotlib.font_manager import FontProperties
from matplotlib import rcParams
from itertools import chain
rc('font', family = 'serif', serif = 'cmr10')
rc('text', usetex=True)
rcParams['text.usetex']=True
rcParams['text.latex.unicode']=True
rcParams['font.family']='serif'
rcParams['font.serif']='cmr10'
station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')
#station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_derived_parameters.npy')
lat_min=-10.
lat_max=60.
lon_min=50.
lon_max=110.
match_bad = re.compile(r'9999')
# Combine date and time lists
dt=[]
x=[]
pdt=[]
dts=[]
#plt_st=np.empty(len(station_data), dtype=str)
#plt_tim=np.zeros((len(station_data)))
plt_tim=[]
plt_st=[]
s= sorted(station_data, key=lambda station: station[0])
#plt_st,plt_tim=zip(*s)
# Convert to python date time
for l,q in enumerate(s):
if float(q[1]) > lat_min and float(q[1]) < lat_max and float(q[2]) > lon_min and float(q[2]) < lon_max:
#print q[3][0][0]
try:
match_station= re.compile(r"%s" % q[3][0][0][0:5])
except Exception,e:
#print e
pass
if match_station.search(str(station_list_pick)) is not None:
#print q[1]
times=[]
plt_st.append(q[0].lower().title())
for n,m in enumerate(q[3]):
if match_bad.search(m[0]) is None:
if match_bad.search(m[1]) is None:
p = datetime.datetime.strptime('%s%s' % (m[0][5:16], m[1]), "%Y%m%d%H%M")
#print p
if p>=date_min and p<=date_max:
times.append(p)
#if match_bad.search(m[1]) is not None:
#p = datetime.datetime.strptime('%s%s' % (m[0][5:16], m[2]), "%Y%m%d%H")
#print p
#if p>=date_min and p<=date_max:
# times.append(p)
plt_tim.append(times)
#print plt_tim
#print len(plt_tim)
#print min(list(chain(*plt_tim)))
b = np.arange(10)
ys = [i+b+(i*b)**2 for i in range(l+1)]
colors = cm.rainbow(np.linspace(0, 1, len(ys)))
time_min = min(list(chain(*plt_tim)))
time_max = max(list(chain(*plt_tim)))
plt.figure(figsize=(18,12))
plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m/%d/%H:%M'))
plt.gca().xaxis.set_major_locator(mdates.DayLocator(interval=4))
plt.gca().xaxis.set_minor_locator(mdates.DayLocator(interval=1))
plt.gca().yaxis.set_ticks([])
#gridlines=plt.gca().get_xgridlines()
plt.gca().xaxis.grid(True)
#for line in gridlines:
# line.set_linestyle('-')
pos_count=0
for t, x in enumerate(plt_tim):
if x!=[]:
y = [1*pos_count] * len(x)
pos_count+= 1
c = colors[t]
l=plt_st[t].replace('_',' ')
plt.scatter(x, y, label=l, color = c)
plt.gcf().autofmt_xdate()
#print x
pos = [min(x), (y[0])]
plt.text(time_min-datetime.timedelta(hours=8), pos[1], l, size=14, rotation=0, ha="right", va="center")
#plt.legend(loc='upper left', prop={'size':4})
plt.ylim([-0.5,y[0]+0.5])
plt.xlim([time_min-datetime.timedelta(hours=5), time_max+datetime.timedelta(hours=5)])
plt.title('Sounding date and time for each station in model domain (Measured)', fontsize=14)
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/igra_freq_picked_stations.png', format='png', bbox_inches='tight')
In [2]:
#############
# Count number of strings in list that are the same
#
#
#############
station_list_pick=[43150, 42867, 43014, 42339, 40990, 40948, 41780, 43003, 42182, 42369, 43353, 43295, 43279, 43369, 42867, 42339, 43128, 43285, 43150, 43185, 43371, 43346, 42492, 43192, 42379, 43014, 42027]
from collections import defaultdict
import os
import numpy as np
import re
#import matplotlib
#matplotlib.use('Agg') # Must be before importing matplotlib.pyplot or pylab!
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
station_data=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/Embrace_Period_India_Station_and_sounding_Info_measured_derived.npy')
#match_bad = re.compile(r'9999')
station_list='/nfs/a90/eepdw/Data/Observations/Radiosonde_downloaded_from_NOAA_GUAN/igra-stations.txt'
st_lat_c=[]
st_lon_c=[]
st_nam_c=[]
st_cnt_c=[]
for i, c in enumerate(station_data):
if c[3]:
match_station= re.compile(r"%s" % c[3][0][0][0:5])
if match_station.search(str(station_list_pick)) is not None:
#print st_lat_c
st_lat_c.append(c[1])
st_lon_c.append(c[2])
st_nam_c.append(c[0].lower().title())
#st_cnt_c[i]=str()
#st_namcnt_c[i]=st_nam_c[i] + ': ' + st_cnt_c[i]
#st_namcnt_c.append(st_nam_c[i])
#print station_data
#PLOT TIME AND DATE vS FREQ ALL STATIONS
#print st_nam_c
#print st_lon_c
# create figure and axes instances
fig = plt.figure(figsize=(16,16))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# create polar stereographic Basemap instance.
#m = Basemap(projection='ste',lon_0=lon_mid,lat_0=lat_mid,lat_ts=lat_mid,\
m = Basemap(projection='cyl',\
llcrnrlat=-10.,urcrnrlat=40.,\
llcrnrlon=60.,urcrnrlon=100.,\
rsphere=6371200.,resolution='l',area_thresh=10000)
# draw coastlines, state and country boundaries, edge of map.
m.drawcoastlines()
# draw parallels.
parallels = np.arange(0.,90,10.)
m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
# draw meridians
meridians = np.arange(0.,360.,10.)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
x, y = m(st_lon_c,st_lat_c)
m.scatter(x,y,3,marker='o',color='red')
for i,j,s in zip(x, y, st_nam_c):
plt.text(i, j, '$%s$' % s, fontsize=14)
plt.title('Position and number of soundings in August and September 2011')
plt.savefig('/nfs/a90/eepdw/Figures/Observation_Plots/sounding_station_map_igra_guan_favourites.png', format='png', bbox_inches='tight')
#plt.show()
for x in sorted(st_nam_c):
print x
In [89]:
# Vizag to Afghanistan
station_list_cs=[43150, 43014, 42867, 42339, 40948, 40990]
#station_list=[42027]
#Indexes - pressures = 1, ptemps = 2, vptemps = 3, rh = 4, uwind = 5, vwind = 6, calculated geopotential heights = 7 . . . may change
# Single level indexes - dates_for_plotting_single = 0, pwater = 1, lclpress = 2, lfcpress = 3, lnbpress = 4, cape = 5,cin = 6
In [105]:
import scipy.interpolate
def interpolate_station(pressures, plot_variable, x_points):
interp = scipy.interpolate.interp1d(pressures, plot_variable, bounds_error=False, fill_value=np.nan)
y_interp = interp(x_points)
return y_interp
In [705]:
#
# Loading individual stations numpy files, sorting into datetime bins
# AND then interpolating for each station
# BEFORE appending so that all 1-D interpolated values for stations are in one array with distances for gridding and plotting
# Plot variable index must be changed in script (plot_variable)
from dateutil.relativedelta import relativedelta
from math import sin, cos, atan2, radians, sqrt
from collections import defaultdict
import collections
import bisect
from scipy.interpolate import spline
from scipy.interpolate import UnivariateSpline
#import scipy.interpolate.interp1d
Cross_section_title = 'Vizag_to_Afghanistan'
first_station = 43150
first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
#print first_station_lon
#print first_station_lat
#delta = relativedelta(months=+1)
delta = relativedelta(days=+1)
#print delta
date_min=datetime.datetime(2011,5,25,0,0,0)
date_max=datetime.datetime(2011,7,10,0,0,0)
print delta
grid=[]
d = date_min
while (d <= date_max):
grid.append(d)
d += delta
#print grid
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
pressuress=[[] for i in range(len(grid)+1)]
distancess=[[] for i in range(len(grid)+1)]
all_stat_date_range_pressure_interp=[[] for i in range(len(grid)+1)]
for stat in station_list_cs:
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
station_title, station_lon, station_lat = station_info_search(stat)
dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
# Pressure Levels
pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s.npy' % (stat, date_title))
for d, dates in enumerate(pressure_file[0]):
idx=bisect.bisect(grid,dates)
bins[idx].append(pressure_file[:,d])
for i in bins:
# if
# if bins[i]:
pv_sort = np.array(bins[i]).T
pv_sort = np.sort(pv_sort, axis=1)
plot_variable = np.array(pv_sort[4], dtype=float)
pressures = np.array(pv_sort[1], dtype=float)
y_points=np.linspace(min(pressures), max(pressures), 1000)
weights = np.ones(plot_variable.shape)
mask = np.isnan(plot_variable)
weights[mask] = 0
nan_mask = np.ma.masked_array(np.array(plot_variable, dtype=float), np.isnan(np.array(plot_variable, dtype=float)))
# print len(pressures)
sm = (len(pressures) * np.var(nan_mask))
#sm=len(weights)
# print np.var(nan_mask)
# print sm
# print i
# print pv_sort
plot_variable[mask] = 0 # arbitrary
if sm:
try:
s = UnivariateSpline(pressures, np.array(plot_variable, dtype=float), s=sm, w=weights)
single_station_single_date_range_pressure_interp = s(y_points)
all_stat_date_range_pressure_interp[i].extend(list(single_station_single_date_range_pressure_interp))
pressuress[i].extend(list(y_points))
distancess[i].extend(list([dist_from_first_station for number in (range(len(y_points)))]))
except:
pass
all_stat_date_range_pressure_interp_np=np.array((distancess,pressuress,all_stat_date_range_pressure_interp))
#for i in bins:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section__IND_INTERP_%s_PRESSURES_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , all_stat_date_range_pressure_interp_np[:,i])
#for i in bins_single:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_%s_SINGLE_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , np.array(bins_single[i]))
In [712]:
#
# Loading individual stations numpy files, sorting into datetime bins
# AND then interpolating for each station
# BEFORE appending so that all 1-D interpolated values for stations are in one array with distances for gridding and plotting
# Plot variable index must be changed in script (plot_variable)
from dateutil.relativedelta import relativedelta
from math import sin, cos, atan2, radians, sqrt
from collections import defaultdict
import collections
import bisect
from scipy.interpolate import spline
#import scipy.interpolate.interp1d
Cross_section_title = 'Vizag_to_Afghanistan'
first_station = 43150
first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
#print first_station_lon
#print first_station_lat
#delta = relativedelta(months=+1)
delta = relativedelta(days=+1)
print delta
all_stat_date_range_pressure_interp=[]
grid=[]
d = date_min
while (d <= date_max):
grid.append(d)
d += delta
#print grid
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
pressuress=[[] for i in range(len(grid)+1)]
distancess=[[] for i in range(len(grid)+1)]
all_stat_date_range_pressure_interp=[[] for i in range(len(grid)+1)]
#all_stat_date_range_pressure_interp_np = np.empty((3, len(grid), len(x_points)), dtype=float)
#all_stat_date_range_pressure_interp_dict={}
for stat in station_list_cs:
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
station_title, station_lon, station_lat = station_info_search(stat)
dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
#print dist_from_first_station
# Pressure Levels
pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s.npy' % (stat, date_title))
for d, dates in enumerate(pressure_file[0]):
idx=bisect.bisect(grid,dates)
bins[idx].append(pressure_file[:,d])
for i in bins:
#bins = np.sort(bins, axis=1)
#print i
pv_sort = np.array(bins[i]).T
pv_sort = np.sort(pv_sort, axis=1)
plot_variable = np.array(pv_sort[4], dtype=float)
nan_mask = np.ma.masked_array(np.array(plot_variable, dtype=float), np.isnan(np.array(plot_variable, dtype=float)))
pressures = np.array(pv_sort[1], dtype=float)
y_points=np.linspace(min(pressures), max(pressures), 25)
single_station_single_date_range_pressure_interp = interpolate_station(pressures, nan_mask, y_points) # interpolating at a station in a date range
all_stat_date_range_pressure_interp[i].extend(list(single_station_single_date_range_pressure_interp))
pressuress[i].extend(list(y_points))
distancess[i].extend(list([dist_from_first_station for number in (range(len(y_points)))]))
all_stat_date_range_pressure_interp1d=np.array((distancess,pressuress,all_stat_date_range_pressure_interp))
#for i in bins:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section__IND_INTERP_%s_PRESSURES_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , all_stat_date_range_pressure_interp_np[:,i])
#for i in bins_single:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_%s_SINGLE_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , np.array(bins_single[i]))
In [847]:
# Loading individual stations numpy files, sorting into datetime bins
# AND then interpolating for each station
# BEFORE appending so that all 1-D interpolated values for stations are in one array with distances for gridding and plotting
# Plot variable index must be changed in script (plot_variable)
#Indexes - pressures = 1, ptemps = 2, vptemps = 3, rh = 4, uwind = 5, vwind = 6, calculated geopotential heights = 7 ,
#temps=8, actual vapour pressure=9, wvmr = 10, theta_e = 11, theta_es = 12. . . may change
variable_list={'pressures': 1, 'theta':2, 'vptemps':3, 'rel_hum':4, 'uwind':5, 'vwind':6, 'calculated geopotential heights':7 ,
'temps':8, 'actual vapour pressure':9, 'wvmr':10, 'theta_e':11, 'theta_es':12}
from dateutil.relativedelta import relativedelta
from math import sin, cos, atan2, radians, sqrt
from collections import defaultdict
import collections
import bisect
from scipy.interpolate import spline
#import scipy.interpolate.interp1d
Cross_section_title = 'Vizag_to_Afghanistan'
first_station = 43150
first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
#print first_station_lon
#print first_station_lat
#delta = relativedelta(months=+1)
delta = relativedelta(days=+1)
print delta
all_stat_date_range_pressure_interp=[]
grid=[]
d = date_min
while (d <= date_max):
grid.append(d)
d += delta
#print grid
for variable in ('rel_hum', 'wvmr', 'theta', 'theta_e', 'theta_es', 'theta'):
for key, value in variable_list.iteritems(): # iter on both keys and values
if key.startswith('%s' % variable):
arr_index_var=value
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
pressuress=[[] for i in range(len(grid)+1)]
distancess=[[] for i in range(len(grid)+1)]
all_stat_date_range_pressure_interp=[[] for i in range(len(grid)+1)]
#all_stat_date_range_pressure_interp_np = np.empty((3, len(grid), len(x_points)), dtype=float)
#all_stat_date_range_pressure_interp_dict={}
for stat in station_list_cs:
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
station_title, station_lon, station_lat = station_info_search(stat)
dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
#print dist_from_first_station
# Pressure Levels
#pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s.npy' % (stat, date_title))
pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s_to_%s.npy'\
% (stat, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d') ))
#print '/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s_to_%s.npy'\
# % (stat, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d') )
#print pressure_file
for d, dates in enumerate(pressure_file[0]):
idx=bisect.bisect(grid,dates)
bins[idx].append(pressure_file[:,d])
for i in bins:
#bins = np.sort(bins, axis=1)
#print i
pv_sort = np.array(bins[i]).T
pv_sort = np.sort(pv_sort, axis=1)
plot_variable = np.array(pv_sort[arr_index_var], dtype=float)
nan_mask = np.ma.masked_array(np.array(plot_variable, dtype=float), np.isnan(np.array(plot_variable, dtype=float)))
pressures = np.array(pv_sort[1], dtype=float)
y_points=np.linspace(min(pressures), max(pressures), 25)
single_station_single_date_range_pressure_interp = interpolate_station(pressures, nan_mask, y_points) # interpolating at a station in a date range
all_stat_date_range_pressure_interp[i].extend(list(single_station_single_date_range_pressure_interp))
pressuress[i].extend(list(y_points))
distancess[i].extend(list([dist_from_first_station for number in (range(len(y_points)))]))
all_stat_date_range_pressure_interp_np=np.array((distancess,pressuress,all_stat_date_range_pressure_interp))
np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_%s_%s_%s_%s' \
% (Cross_section_title, variable, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta)\
, all_stat_date_range_pressure_interp_np)
#for i in bins:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section__IND_INTERP_%s_PRESSURES_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , all_stat_date_range_pressure_interp_np[:,i])
#for i in bins_single:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_%s_SINGLE_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , np.array(bins_single[i]))
In [847]:
# Loading individual stations numpy files, sorting into datetime bins
# AND then interpolating for each station
# BEFORE appending so that all 1-D interpolated values for stations are in one array with distances for gridding and plotting
# Plot variable index must be changed in script (plot_variable)
#Indexes - pressures = 1, ptemps = 2, vptemps = 3, rh = 4, uwind = 5, vwind = 6, calculated geopotential heights = 7 ,
#temps=8, actual vapour pressure=9, wvmr = 10, theta_e = 11, theta_es = 12. . . may change
variable_list={'pressures': 1, 'theta':2, 'vptemps':3, 'rel_hum':4, 'uwind':5, 'vwind':6, 'calculated geopotential heights':7 ,
'temps':8, 'actual vapour pressure':9, 'wvmr':10, 'theta_e':11, 'theta_es':12}
from dateutil.relativedelta import relativedelta
from math import sin, cos, atan2, radians, sqrt
from collections import defaultdict
import collections
import bisect
from scipy.interpolate import spline
#import scipy.interpolate.interp1d
Cross_section_title = 'Vizag_to_Afghanistan'
first_station = 43150
first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
#print first_station_lon
#print first_station_lat
#delta = relativedelta(months=+1)
delta = relativedelta(days=+1)
print delta
all_stat_date_range_pressure_interp=[]
grid=[]
d = date_min
while (d <= date_max):
grid.append(d)
d += delta
#print grid
for variable in ('rel_hum', 'wvmr', 'theta', 'theta_e', 'theta_es', 'theta'):
for key, value in variable_list.iteritems(): # iter on both keys and values
if key.startswith('%s' % variable):
arr_index_var=value
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
pressuress=[[] for i in range(len(grid)+1)]
distancess=[[] for i in range(len(grid)+1)]
all_stat_date_range_pressure_interp=[[] for i in range(len(grid)+1)]
#all_stat_date_range_pressure_interp_np = np.empty((3, len(grid), len(x_points)), dtype=float)
#all_stat_date_range_pressure_interp_dict={}
for stat in station_list_cs:
bins=collections.defaultdict(list)
bins_single=collections.defaultdict(list)
station_title, station_lon, station_lat = station_info_search(stat)
dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
#print dist_from_first_station
# Pressure Levels
#pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s.npy' % (stat, date_title))
pressure_file = np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s_to_%s.npy'\
% (stat, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d') ))
#print '/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Single_Station_PRESSURES_%s_%s_to_%s.npy'\
# % (stat, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d') )
#print pressure_file
for d, dates in enumerate(pressure_file[0]):
idx=bisect.bisect(grid,dates)
bins[idx].append(pressure_file[:,d])
for i in bins:
#bins = np.sort(bins, axis=1)
#print i
pv_sort = np.array(bins[i]).T
pv_sort = np.sort(pv_sort, axis=1)
plot_variable = np.array(pv_sort[arr_index_var], dtype=float)
nan_mask = np.ma.masked_array(np.array(plot_variable, dtype=float), np.isnan(np.array(plot_variable, dtype=float)))
pressures = np.array(pv_sort[1], dtype=float)
y_points=np.linspace(min(pressures), max(pressures), 25)
single_station_single_date_range_pressure_interp = interpolate_station(pressures, nan_mask, y_points) # interpolating at a station in a date range
all_stat_date_range_pressure_interp[i].extend(list(single_station_single_date_range_pressure_interp))
pressuress[i].extend(list(y_points))
distancess[i].extend(list([dist_from_first_station for number in (range(len(y_points)))]))
all_stat_date_range_pressure_interp_np=np.array((distancess,pressuress,all_stat_date_range_pressure_interp))
np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_%s_%s_%s_%s' \
% (Cross_section_title, variable, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta)\
, all_stat_date_range_pressure_interp_np)
#for i in bins:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section__IND_INTERP_%s_PRESSURES_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , all_stat_date_range_pressure_interp_np[:,i])
#for i in bins_single:
# np.save('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/Radiosonde_Cross_Section_%s_SINGLE_%s' \
# % (Cross_section_title, grid[i].strftime('%Y%m'))\
# , np.array(bins_single[i]))
In [752]:
#Indexes - pressures = 1, ptemps = 2, vptemps = 3, rh = 4, uwind = 5, vwind = 6, calculated geopotential heights = 7 , temps=8, actual vapour pressure=9, wvmr = 10, theta_e = 11, theta_es = 12. . . may change
# Single level indexes - dates_for_plotting_single = 0, pwater = 1, lclpress = 2, lfcpress = 3, lnbpress = 4, cape = 5,cin = 6
def equiv_pot_temp(T_s, theta, r, c_pa):
L = 3136.17 - 2.34 * T_s
theta_e = theta * (1 + r * L / (c_pa * T_s))
return theta_e
def wvmr(act_vap_press, pressure):
WVMR=0.622*(act_vap_press)/((pressure)-(act_vap_press))
return WVMR
#c_pa=1004 # Specific heat of dry air at constant pressure - J/kg/K
#equiv_pot_temp(ps[8]/10, ps[2]/10, WVMR, c_pa)
In [612]:
def grid_data_cs(pressure, distance, param):
xi=np.linspace(0, max(distance), 1000)
yi=np.linspace(min(pressure), max(pressure), 500)
#yi=np.array([1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20,10], dtype=float)
#yi=np.array([10, 20, 30, 50, 70, 100, 150, 200, 250, 300, 400, 500, 700, 850, 925, 1000]*100, dtype=float)
zi = ml.griddata(distance, pressure,param,xi, yi, interp='nn')
#zi = scipy.interpolate.griddata((matplotlib.dates.date2num(q[0]),q[1]), param, (xi[None,:],yi[:,None]), method='linear')
return xi,yi,zi
In [590]:
def station_name_plot(station_list_cs, first_station):
y_offset_text=0
first_station_title, first_station_lon, first_station_lat = station_info_search(first_station)
for stat in station_list_cs:
station_title, station_lon, station_lat = station_info_search(stat)
dist_from_first_station = calculate_distance_from_first_station(stat, first_station_lon, first_station_lat, station_lat, station_lon)
plt.axvline(x=dist_from_first_station, ymin=0, ymax=1, label=station_title, color='k')
plt.text(dist_from_first_station+0.1,max(yi)/100+20,station_title,rotation=-45)
y_offset_text=+1
In [706]:
first_station=43150
for i in bins:
nan_mask = np.ma.masked_array(np.array(all_stat_date_range_pressure_interp_np[2,i], dtype=float), np.isnan(np.array(all_stat_date_range_pressure_interp_np[2,i], dtype=float)))
xi,yi,zi = grid_data_cs(np.array(all_stat_date_range_pressure_interp_np[1,i], dtype=float), np.array(all_stat_date_range_pressure_interp_np[0,i], dtype=float), nan_mask/10)
#print zi
cont,cbar = plot_rad_cs(xi,yi,zi)
station_name_plot(station_list_cs, first_station)
cbar.set_label('\%', rotation=90)
plt.title('%s %s Cross-Section of Relative Humidity from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
plt.show()
In [707]:
def grid_data_cs(pressure, distance, param):
xi=np.linspace(0, max(distance), 1000)
yi=np.linspace(min(pressure), max(pressure), 25)
#yi=np.array([1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100, 70, 50, 30, 20,10], dtype=float)
#yi=np.array([10, 20, 30, 50, 70, 100, 150, 200, 250, 300, 400, 500, 700, 850, 925, 1000]*100, dtype=float)
zi = ml.griddata(distance, pressure,param,xi, yi, interp='nn')
#zi = scipy.interpolate.griddata((matplotlib.dates.date2num(q[0]),q[1]), param, (xi[None,:],yi[:,None]), method='linear')
return xi,yi,zi
In [879]:
#Daily
# Remember to change savefig filename depending on timedelta - daily, weekly, monthly
first_station=43150
Cross_section_title='Vizag_to_Afghanistan'
datetime_string_format='%Y%m%d'
# Relative Humidity
all_stat_date_range_pressure_interp1d=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_rel_hum_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%s' % datetime_string_format), date_max.strftime('%Y%m%d'), delta))
for i in bins:
ar=all_stat_date_range_pressure_interp1d[:,i]
#arr=np.array(ar, dtype=float)
nan_mask = np.ma.masked_array(ar[2], np.isnan(ar[2]))
xi,yi,zi = grid_data_cs(np.array(ar[1], dtype=float), np.array(ar[0], dtype=float), nan_mask/10)
#print zi
max_contour=100
min_contour=0
cont,cbar = plot_rad_cs(xi,yi,zi, min_contour, max_contour)
station_name_plot(station_list_cs, first_station)
cbar.set_label('\%', rotation=90)
plt.title('%s %s Cross-Section of Relative Humidity from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Relative_Humidity.png' % (C_section, year, bins[i][0][0].strftime("%Y%m%d")), format='png', bbox_inches='tight')
#Theta E
all_stat_date_range_pressure_interp1d=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_theta_e_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta))
for i in bins:
ar=all_stat_date_range_pressure_interp1d[:,i]
#arr=np.array(ar, dtype=float)
nan_mask = np.ma.masked_array(ar[2], np.isnan(ar[2]))
xi,yi,zi = grid_data_cs(np.array(ar[1], dtype=float), np.array(ar[0], dtype=float), nan_mask)
#print zi
max_contour=400
min_contour=300
cont,cbar = plot_rad_cs(xi,yi,zi, min_contour, max_contour)
station_name_plot(station_list_cs, first_station)
cbar.set_label('K', rotation=90)
plt.title('%s %s Cross-Section of Equivalent Potential Temperature from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Theta_E.png' % (C_section, year, bins[i][0][0].strftime("%Y%m%d")), format='png', bbox_inches='tight')
#Theta Es
all_stat_date_range_pressure_interp1d=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_theta_es_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta))
for i in bins:
ar=all_stat_date_range_pressure_interp1d[:,i]
#arr=np.array(ar, dtype=float)
nan_mask = np.ma.masked_array(ar[2], np.isnan(ar[2]))
xi,yi,zi = grid_data_cs(np.array(ar[1], dtype=float), np.array(ar[0], dtype=float), nan_mask)
#print zi
max_contour=400
min_contour=300
cont,cbar = plot_rad_cs(xi,yi,zi, min_contour, max_contour)
station_name_plot(station_list_cs, first_station)
cbar.set_label('K', rotation=90)
plt.title('%s %s Cross-Section of Saturated Equivalent Potential Temperature from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Theta_E_S.png' % (C_section, year, bins[i][0][0].strftime("%Y%m%d")), format='png', bbox_inches='tight')
# WVMR
all_stat_date_range_pressure_interp1d=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_wvmr_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta))
for i in bins:
ar=all_stat_date_range_pressure_interp1d[:,i]
#arr=np.array(ar, dtype=float)
nan_mask = np.ma.masked_array(ar[2], np.isnan(ar[2]))
xi,yi,zi = grid_data_cs(np.array(ar[1], dtype=float), np.array(ar[0], dtype=float), nan_mask*1000)
#print zi
max_contour=15
min_contour=0
cont,cbar = plot_rad_cs(xi,yi,zi, min_contour, max_contour)
station_name_plot(station_list_cs, first_station)
cbar.set_label('g/kg', rotation=90)
plt.title('%s %s Cross-Section of Water Vapour Mixing Ratio from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
#plt.show()
plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_WVMR.png' % (C_section, year, bins[i][0][0].strftime("%Y%m%d")), format='png', bbox_inches='tight')
#Theta Es - Theta E - Buoyancy
theta_e=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_theta_e_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta))
theta_es=np.load('/nfs/a90/eepdw/Data/Observations/Radiosonde_Numpy/'\
'Radiosonde_Cross_Section_IND_INTERP_%s_PRESSURES_theta_es_%s_%s_%s.npy' \
% (Cross_section_title, date_min.strftime('%Y%m%d'), date_max.strftime('%Y%m%d'), delta))
all_stat_date_range_pressure_interp1d=theta_e-theta_es
for i in bins:
ar=all_stat_date_range_pressure_interp1d[:,i]
#arr=np.array(ar, dtype=float)
nan_mask = np.ma.masked_array(ar[2], np.isnan(ar[2]))
xi,yi,zi = grid_data_cs(np.array(ar[1], dtype=float), np.array(ar[0], dtype=float), nan_mask)
#print zi
max_contour=20
min_contour=0
cont,cbar = plot_rad_cs(xi,yi,zi, min_contour, max_contour)
station_name_plot(station_list_cs, first_station)
cbar.set_label('K', rotation=90)
plt.title('%s %s Cross-Section of Saturated Equivalent Potential Temperature from Radiosonde Soundings' % (bins[i][0][0].strftime("%d %B"), C_section.replace('_',' ') ))
#plt.show()
#plt.savefig('/nfs/a90/eepdw/Figures/Radiosonde/Cross_Sections/%s_%s_%s_Theta_E_S_minus_Theta_E.png' % (C_section, year, bins[i][0][0].strftime("%Y%m%d")), format='png', bbox_inches='tight')